Daily_Weather

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(tidyr)
library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(ggplot2)
library(forecast)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(fable)

## Loading required package: fabletools

data <- read.csv("daily_weather.csv")
summary(data)

##      number       air_pressure_9am  air_temp_9am   avg_wind_direction_9am
##  Min.   :   0.0   Min.   :908.0    Min.   :36.75   Min.   : 15.50        
##  1st Qu.: 273.5   1st Qu.:916.5    1st Qu.:57.28   1st Qu.: 65.97        
##  Median : 547.0   Median :918.9    Median :65.72   Median :166.00        
##  Mean   : 547.0   Mean   :918.9    Mean   :64.93   Mean   :142.24        
##  3rd Qu.: 820.5   3rd Qu.:921.2    3rd Qu.:73.45   3rd Qu.:191.00        
##  Max.   :1094.0   Max.   :929.3    Max.   :98.91   Max.   :343.40        
##                   NA's   :3        NA's   :5       NA's   :4             
##  avg_wind_speed_9am max_wind_direction_9am max_wind_speed_9am
##  Min.   : 0.6935    Min.   : 28.90         Min.   : 1.186    
##  1st Qu.: 2.2488    1st Qu.: 76.55         1st Qu.: 3.067    
##  Median : 3.8713    Median :177.30         Median : 4.944    
##  Mean   : 5.5083    Mean   :148.95         Mean   : 7.020    
##  3rd Qu.: 7.3372    3rd Qu.:201.23         3rd Qu.: 8.948    
##  Max.   :23.5550    Max.   :312.20         Max.   :29.841    
##  NA's   :3          NA's   :3              NA's   :4         
##  rain_accumulation_9am rain_duration_9am relative_humidity_9am
##  Min.   : 0.0000       Min.   :    0.0   Min.   : 6.09        
##  1st Qu.: 0.0000       1st Qu.:    0.0   1st Qu.:15.09        
##  Median : 0.0000       Median :    0.0   Median :23.18        
##  Mean   : 0.2031       Mean   :  294.1   Mean   :34.24        
##  3rd Qu.: 0.0000       3rd Qu.:    0.0   3rd Qu.:45.40        
##  Max.   :24.0200       Max.   :17704.0   Max.   :92.62        
##  NA's   :6             NA's   :3                              
##  relative_humidity_3pm
##  Min.   : 5.30        
##  1st Qu.:17.39        
##  Median :24.38        
##  Mean   :35.34        
##  3rd Qu.:52.06        
##  Max.   :92.25        
## 

data <- na.omit(data)

air_pressure_9am <- data$air_pressure_9am
air_temp_9am <- data$air_temp_9am
avg_wind_direction_9am <- data$avg_wind_direction_9am
avg_wind_speed_9am <- data$avg_wind_speed_9am
max_wind_direction_9am <- data$max_wind_direction_9am
max_wind_speed_9am <- data$max_wind_speed_9am
rain_accumulation_9am <- data$rain_accumulation_9am
rain_duration_9am <- data$rain_duration_9am
relative_humidity_9am <- data$relative_humidity_9am
relative_humidity_3pm <- data$relative_humidity_3pm

N <- length(relative_humidity_3pm)

#Standardizing the data and converting to time series object
air_pressure_9am <- air_pressure_9am - mean(air_pressure_9am)
air_pressure_9am <- ts(air_pressure_9am, start = 1, end = N)

air_temp_9am <- air_temp_9am - mean(air_temp_9am)
air_temp_9am <- ts(air_temp_9am, start = 1, end = N)

avg_wind_direction_9am <- avg_wind_direction_9am - mean(avg_wind_direction_9am)
avg_wind_direction_9am <- ts(avg_wind_direction_9am, start = 1, end = N)

avg_wind_speed_9am <- avg_wind_speed_9am - mean(avg_wind_speed_9am)
avg_wind_speed_9am <- ts(avg_wind_speed_9am, start = 1, end = N)

max_wind_direction_9am <- max_wind_direction_9am - mean(max_wind_direction_9am)
max_wind_direction_9am <- ts(max_wind_direction_9am, start = 1, end = N)

max_wind_speed_9am <- max_wind_speed_9am - mean(max_wind_speed_9am)
max_wind_speed_9am <- ts(max_wind_speed_9am, start = 1, end = N)

rain_accumulation_9am <- rain_accumulation_9am - mean(rain_accumulation_9am)
rain_accumulation_9am <- ts(rain_accumulation_9am, start = 1, end = N)

rain_duration_9am <- rain_duration_9am - mean(rain_duration_9am)
rain_duration_9am <- ts(rain_duration_9am, start = 1, end = N)

relative_humidity_9am <- relative_humidity_9am - mean(relative_humidity_9am)
relative_humidity_9am <- ts(relative_humidity_9am, start = 1, end = N)

relative_humidity_3pm <- relative_humidity_3pm - mean(relative_humidity_3pm)
relative_humidity_3pm <- ts(relative_humidity_3pm, start = 1, end = N)


plot(air_pressure_9am, xlab="days", ylab = "air_pressure", main = "Air Pressure at 9AM")

plot(air_temp_9am, xlab="days", ylab = "air_temperature", main = "Air Temperature at 9AM")

plot(avg_wind_direction_9am, xlab="days", ylab = "avg_wind_direction_9am", main = "Average Wind Direction at 9AM")

plot(avg_wind_speed_9am, xlab="days", ylab = "avg_wind_speed_9am", main = "Average Wind Speed at 9AM")

plot(max_wind_direction_9am, xlab="days", ylab = "max_wind_direction_9am", main = "Max wind direction at 9AM")

plot(max_wind_speed_9am, xlab="days", ylab = "max_wind_speed_9am", main = "Max Wind Speed at 9AM")

plot(rain_accumulation_9am, xlab="days", ylab = "rain_accumulation_9am", main = "Rain Accumulation at 9AM")

plot(rain_duration_9am, xlab="days", ylab = "rain_duration_9am", main = "Rain Duration at 9AM")

plot(relative_humidity_9am, xlab="days", ylab = "relative_humidity_9am", main = "Relative Humidity at 9AM")

plot(relative_humidity_3pm, xlab="days", ylab = "relative_humidity_3pm", main = "Relative Humidity at 3PM")

acf(air_pressure_9am)

acf(air_temp_9am)

acf(avg_wind_direction_9am)

acf(avg_wind_speed_9am)

acf(max_wind_direction_9am)

acf(max_wind_speed_9am)

acf(rain_accumulation_9am)

acf(rain_duration_9am)

acf(relative_humidity_9am)

acf(relative_humidity_3pm)

#Differencing the variables to make it stationary
diff_air_pressure_9am <- diff(air_pressure_9am)
plot(diff_air_pressure_9am, xlab="days", ylab = "air_pressure_9am", main = "Differenced Air Pressure at 9am Plot")

diff_air_temp_9am <- diff(air_temp_9am)
plot(diff_air_temp_9am, xlab="days", ylab = "air_temp_9am", main = "Differenced Air Temperature at 9am Plot")

diff_avg_wind_direction_9am <- diff(avg_wind_direction_9am)
plot(diff_avg_wind_direction_9am, xlab="days", ylab = "avg_wind_direction_9am", main = "Differenced Average Wind Direction at 9am Plot")

diff_avg_wind_speed_9am <- diff(avg_wind_speed_9am)
plot(diff_avg_wind_speed_9am, xlab="days", ylab = "avg_wind_speed_9am", main = "Differenced Average Wind Speed at 9am Plot")

diff_max_wind_direction_9am <- diff(max_wind_direction_9am)
plot(diff_max_wind_direction_9am, xlab="days", ylab = "max_wind_direction_9am", main = "Differenced Mean Reading at Useless Load Plot")

diff_max_wind_speed_9am <- diff(max_wind_speed_9am)
plot(diff_max_wind_speed_9am, xlab="days", ylab = "max_wind_speed_9am", main = "Differenced Max Wind Speed at 9am Plot")

diff_rain_accumulation_9am <- diff(rain_accumulation_9am)
plot(diff_rain_accumulation_9am, xlab="days", ylab = "rain_accumulation_9am", main = "Differenced Rain Accumulation at 9am Plot")

diff_rain_duration_9am <- diff(rain_duration_9am)
plot(diff_rain_duration_9am, xlab="days", ylab = "rain_duration_9am", main = "Differenced Rain Duration at 9am Plot")

diff_relative_humidity_9am <- diff(relative_humidity_9am)
plot(diff_relative_humidity_9am, xlab="days", ylab = "relative_humidity_9am", main = "Differenced Relative Humidity at 9am Plot")

diff_relative_humidity_3pm <- diff(relative_humidity_3pm)
plot(diff_relative_humidity_3pm, xlab="days", ylab = "relative_humidity_3pm", main = "Differenced Relative Humidity at 3pm Plot")

adf.test(diff_air_pressure_9am)

## Warning in adf.test(diff_air_pressure_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_air_pressure_9am
## Dickey-Fuller = -16.148, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_air_temp_9am)

## Warning in adf.test(diff_air_temp_9am): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_air_temp_9am
## Dickey-Fuller = -15.332, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_avg_wind_direction_9am)

## Warning in adf.test(diff_avg_wind_direction_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_avg_wind_direction_9am
## Dickey-Fuller = -17.135, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_avg_wind_speed_9am)

## Warning in adf.test(diff_avg_wind_speed_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_avg_wind_speed_9am
## Dickey-Fuller = -15.752, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_max_wind_direction_9am)

## Warning in adf.test(diff_max_wind_direction_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_max_wind_direction_9am
## Dickey-Fuller = -16.261, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_max_wind_speed_9am)

## Warning in adf.test(diff_max_wind_speed_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_max_wind_speed_9am
## Dickey-Fuller = -15.698, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_rain_accumulation_9am)

## Warning in adf.test(diff_rain_accumulation_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_rain_accumulation_9am
## Dickey-Fuller = -16.832, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_rain_duration_9am)

## Warning in adf.test(diff_rain_duration_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_rain_duration_9am
## Dickey-Fuller = -16.862, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_relative_humidity_9am)

## Warning in adf.test(diff_relative_humidity_9am): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_relative_humidity_9am
## Dickey-Fuller = -16.118, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

adf.test(diff_relative_humidity_3pm)

## Warning in adf.test(diff_relative_humidity_3pm): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_relative_humidity_3pm
## Dickey-Fuller = -15.81, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

#ACF
acf(diff_air_pressure_9am)

acf(diff_air_temp_9am)

acf(diff_avg_wind_direction_9am)

acf(diff_avg_wind_speed_9am)

acf(diff_max_wind_direction_9am)

acf(diff_max_wind_speed_9am)

acf(diff_rain_accumulation_9am)

acf(diff_rain_duration_9am)

acf(diff_relative_humidity_9am)

acf(diff_relative_humidity_3pm)

#PACF
pacf(diff_air_pressure_9am)

pacf(diff_air_temp_9am)

pacf(diff_avg_wind_direction_9am)

pacf(diff_avg_wind_speed_9am)

pacf(diff_max_wind_direction_9am)

pacf(diff_max_wind_speed_9am)

pacf(diff_rain_accumulation_9am)

pacf(diff_rain_duration_9am)

pacf(diff_relative_humidity_9am)

pacf(diff_relative_humidity_3pm)

#CCF
ccf(diff_relative_humidity_3pm, diff_air_pressure_9am)

ccf(diff_relative_humidity_3pm, diff_air_temp_9am)

ccf(diff_relative_humidity_3pm, diff_avg_wind_direction_9am)

ccf(diff_relative_humidity_3pm, diff_avg_wind_speed_9am)

ccf(diff_relative_humidity_3pm, diff_max_wind_direction_9am)

ccf(diff_relative_humidity_3pm, diff_max_wind_speed_9am)

ccf(diff_relative_humidity_3pm, diff_rain_accumulation_9am)

ccf(diff_relative_humidity_3pm, diff_rain_duration_9am)

ccf(diff_relative_humidity_3pm, diff_relative_humidity_9am)

spectrum(diff_air_pressure_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Air Pressure at 9am')

spectrum(diff_air_temp_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Air Temperature at 9am')

spectrum(diff_avg_wind_direction_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Average Wind Direction at 9am')

spectrum(diff_avg_wind_speed_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Average Wind Speed at 9am')

spectrum(diff_max_wind_direction_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Max Wind Direction at 9am')

spectrum(diff_max_wind_speed_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Max Wind Speed at 9am')

spectrum(diff_rain_accumulation_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Rain Accumulation at 9am')

spectrum(diff_rain_duration_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Rain Duration at 9am')

spectrum(diff_relative_humidity_9am, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Relative Humidity at 9am')

spectrum(diff_relative_humidity_3pm, log='no', xlab='Frequencies',
         ylab='Power',
         main='Periodogram Relative Humidity at 3pm')

suppressWarnings({

#Based on observations from ACF and CCF creating new lagged inputs and output
air_pressure_9am_lag1 <- filter(diff_air_pressure_9am, c(0,1), sides=1)
air_pressure_9am_lag2 <- filter(air_pressure_9am_lag1, c(0,1), sides=1)

air_temp_9am_lag1 <- filter(diff_air_temp_9am, c(0,1), sides=1)
air_temp_9am_lag2 <- filter(air_temp_9am_lag1, c(0,1), sides=1)

avg_wind_direction_9am_lag1 <- filter(diff_avg_wind_direction_9am, c(0,1), sides=1)
avg_wind_direction_9am_lag2 <- filter(avg_wind_direction_9am_lag1, c(0,1), sides=1)

avg_wind_speed_9am_lag1 <- filter(diff_avg_wind_speed_9am, c(0,1), sides=1)
avg_wind_speed_9am_lag2 <- filter(avg_wind_speed_9am_lag1, c(0,1), sides=1)

max_wind_direction_9am_lag1 <- filter(diff_max_wind_direction_9am, c(0,1), sides=1)
max_wind_direction_9am_lag2 <- filter(max_wind_direction_9am_lag1, c(0,1), sides=1)

max_wind_speed_9am_lag1 <- filter(diff_max_wind_speed_9am, c(0,1), sides=1)
max_wind_speed_9am_lag2 <- filter(max_wind_speed_9am_lag1, c(0,1), sides=1)

rain_accumulation_9am_lag1 <- filter(diff_rain_accumulation_9am, c(0,1), sides=1)
rain_accumulation_9am_lag2 <- filter(rain_accumulation_9am_lag1, c(0,1), sides=1)

rain_duration_9am_lag1 <- filter(diff_rain_duration_9am, c(0,1), sides=1)
rain_duration_9am_lag2 <- filter(rain_duration_9am_lag1, c(0,1), sides=1)

relative_humidity_9am_lag1 <- filter(diff_relative_humidity_9am, c(0,1), sides=1)
relative_humidity_9am_lag2 <- filter(relative_humidity_9am_lag1, c(0,1), sides=1)

relative_humidity_3pm_lag1 <- filter(diff_relative_humidity_3pm, c(0,1), sides=1)
relative_humidity_3pm_lag2 <- filter(relative_humidity_3pm_lag1, c(0,1), sides=1)

diff_air_pressure_9am <- diff_air_pressure_9am[1:N]
air_pressure_9am_lag1 <- air_pressure_9am_lag1[1:N]
air_pressure_9am_lag2 <- air_pressure_9am_lag2[1:N]

diff_air_temp_9am <- diff_air_temp_9am[1:N]
air_temp_9am_lag1 <- air_temp_9am_lag1[1:N]
air_temp_9am_lag2 <- air_temp_9am_lag2[1:N]

diff_avg_wind_direction_9am <- diff_avg_wind_direction_9am[1:N]
avg_wind_direction_9am_lag1 <- avg_wind_direction_9am_lag1[1:N]
avg_wind_direction_9am_lag2 <- avg_wind_direction_9am_lag2[1:N]

diff_avg_wind_speed_9am <- diff_avg_wind_speed_9am[1:N]
avg_wind_speed_9am_lag1 <- avg_wind_speed_9am_lag1[1:N]
avg_wind_speed_9am_lag2 <- avg_wind_speed_9am_lag2[1:N]

diff_max_wind_direction_9am <- diff_max_wind_direction_9am[1:N]
max_wind_direction_9am_lag1 <- max_wind_direction_9am_lag1[1:N]
max_wind_direction_9am_lag2 <- max_wind_direction_9am_lag2[1:N]

diff_max_wind_speed_9am <- diff_max_wind_speed_9am[1:N]
max_wind_speed_9am_lag1 <- max_wind_speed_9am_lag1[1:N]
max_wind_speed_9am_lag2 <- max_wind_speed_9am_lag2[1:N]

diff_rain_accumulation_9am <- diff_rain_accumulation_9am[1:N]
rain_accumulation_9am_lag1 <- rain_accumulation_9am_lag1[1:N]
rain_accumulation_9am_lag2 <- rain_accumulation_9am_lag2[1:N]

diff_rain_duration_9am <- diff_rain_duration_9am[1:N]
rain_duration_9am_lag1 <- rain_duration_9am_lag1[1:N]
rain_duration_9am_lag2 <- rain_duration_9am_lag2[1:N]

diff_relative_humidity_9am <- diff_relative_humidity_9am[1:N]
relative_humidity_9am_lag1 <- relative_humidity_9am_lag1[1:N]
relative_humidity_9am_lag2 <- relative_humidity_9am_lag2[1:N]

diff_relative_humidity_3pm <- diff_relative_humidity_3pm[1:N]
relative_humidity_3pm_lag1 <- relative_humidity_3pm_lag1[1:N]
relative_humidity_3pm_lag2 <- relative_humidity_3pm_lag2[1:N]

time_prd <- 1:N-1
sin_pred <- sin(2 * pi * time_prd * (1/365))
cos_pred <- cos(2 * pi * time_prd * (1/365))

data_arima2 <- cbind(air_pressure_9am_lag1, air_pressure_9am_lag2,
                     air_temp_9am_lag1, air_temp_9am_lag2, 
                     avg_wind_direction_9am_lag1, avg_wind_direction_9am_lag2,
                     avg_wind_speed_9am_lag1, avg_wind_speed_9am_lag2, 
                     rain_duration_9am_lag1, rain_duration_9am_lag2,
                     relative_humidity_9am_lag1, relative_humidity_9am_lag2,
                     relative_humidity_3pm_lag1, relative_humidity_3pm_lag2, sin_pred, cos_pred)

data_arima1 <- cbind(air_pressure_9am_lag1, air_temp_9am_lag1, avg_wind_direction_9am_lag1, avg_wind_speed_9am_lag1, rain_duration_9am_lag1, relative_humidity_9am_lag1, relative_humidity_3pm_lag1, sin_pred, cos_pred)

})

for(i in 1:5) {
  for (j in 1:5){
    complex_model_lag2 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,], order=c(i,0,j))
    pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])
    
    ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"))
    legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
    
    cat()
    cat("AIC:", AIC(complex_model_lag2), "BIC:", BIC(complex_model_lag2), "p:", i, "q:", j,"\n")
    cat()
  }
}

## AIC: 9065.332 BIC: 9163.447 p: 1 q: 1

## AIC: 9065.844 BIC: 9168.865 p: 1 q: 2

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9064.37 BIC: 9172.296 p: 1 q: 3

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9067.453 BIC: 9180.285 p: 1 q: 4

## AIC: 9070.327 BIC: 9188.065 p: 1 q: 5

## AIC: 9066.041 BIC: 9169.062 p: 2 q: 1

## AIC: 9067.595 BIC: 9175.521 p: 2 q: 2

## AIC: 9066.132 BIC: 9178.964 p: 2 q: 3

## AIC: 9065.692 BIC: 9183.431 p: 2 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9066.021 BIC: 9188.665 p: 2 q: 5

## AIC: 9067.933 BIC: 9175.86 p: 3 q: 1

## Warning in log(s2): NaNs produced

## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced

## AIC: 9068.679 BIC: 9181.512 p: 3 q: 2

## AIC: 9065.597 BIC: 9183.335 p: 3 q: 3

## AIC: 9067.818 BIC: 9190.462 p: 3 q: 4

## AIC: 9068.575 BIC: 9196.124 p: 3 q: 5

## AIC: 9069.878 BIC: 9182.711 p: 4 q: 1

## Warning in log(s2): NaNs produced

## AIC: 9070.332 BIC: 9188.07 p: 4 q: 2

## AIC: 9068.148 BIC: 9190.792 p: 4 q: 3

## AIC: 9069.889 BIC: 9197.439 p: 4 q: 4

## AIC: 9060.747 BIC: 9193.203 p: 4 q: 5

## AIC: 9071.647 BIC: 9189.385 p: 5 q: 1

## AIC: 9071.674 BIC: 9194.318 p: 5 q: 2

## AIC: 9069.755 BIC: 9197.304 p: 5 q: 3

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9064.824 BIC: 9197.279 p: 5 q: 4

## AIC: 9069.728 BIC: 9207.09 p: 5 q: 5

for(i in 1:5) {
  for (j in 1:5){
                        # Fit ARMAX model
                        armax_model <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(i, 0, j), include.mean = FALSE)

                        # Predict using ARMAX model
                        pr_armax <- predict(armax_model, newxreg = data_arima1[1001:1050,])

                        # Plot actual vs predicted values
                        ts.plot(diff_relative_humidity_3pm[1000:1049], pr_armax$pred, lty=1:2, col=c("blue", "red"))
                        legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
                        cat()
                        # Print AIC and BIC scores
                        cat("AIC:", AIC(armax_model), "BIC:", BIC(armax_model), "p:", i, "q:", j,"\n")
                        cat()
  }
}

## AIC: 9069.574 BIC: 9128.455 p: 1 q: 1

## AIC: 9071.266 BIC: 9135.054 p: 1 q: 2

## AIC: 9073.089 BIC: 9141.783 p: 1 q: 3

## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced

## AIC: 9074.745 BIC: 9148.346 p: 1 q: 4

## AIC: 9075.388 BIC: 9153.896 p: 1 q: 5

## AIC: 9071.029 BIC: 9134.817 p: 2 q: 1

## AIC: 9073.199 BIC: 9141.894 p: 2 q: 2

## AIC: 9074.595 BIC: 9148.196 p: 2 q: 3

## AIC: 9070.652 BIC: 9149.16 p: 2 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9072.456 BIC: 9155.87 p: 2 q: 5

## AIC: 9072.848 BIC: 9141.542 p: 3 q: 1

## AIC: 9073.779 BIC: 9147.38 p: 3 q: 2

## AIC: 9075.753 BIC: 9154.261 p: 3 q: 3

## AIC: 9072.472 BIC: 9155.887 p: 3 q: 4

## AIC: 9072.756 BIC: 9161.078 p: 3 q: 5

## AIC: 9074.385 BIC: 9147.986 p: 4 q: 1

## AIC: 9075.912 BIC: 9154.42 p: 4 q: 2

## AIC: 9077.974 BIC: 9161.389 p: 4 q: 3

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9073.77 BIC: 9162.092 p: 4 q: 4

## AIC: 9073.691 BIC: 9166.919 p: 4 q: 5

## AIC: 9075.776 BIC: 9154.284 p: 5 q: 1

## AIC: 9077.228 BIC: 9160.643 p: 5 q: 2

## AIC: 9074.709 BIC: 9163.03 p: 5 q: 3

## AIC: 9076.222 BIC: 9169.45 p: 5 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9069.652 BIC: 9167.787 p: 5 q: 5

for(i in 1:5) {
  for (j in 1:5){
    complex_model_lag1 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(i,0,j))
    pr <- predict(complex_model_lag1, newxreg = data_arima1[1001:1050,])
    
    ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"))
    legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
    
    cat()
    cat("AIC:", AIC(complex_model_lag1), "BIC:", BIC(complex_model_lag1), "p:", i, "q:", j,"\n")
    cat()
  }
}

## AIC: 9070.219 BIC: 9134.007 p: 1 q: 1

## AIC: 9071.847 BIC: 9140.541 p: 1 q: 2

## AIC: 9073.646 BIC: 9147.247 p: 1 q: 3

## AIC: 9075.057 BIC: 9153.566 p: 1 q: 4

## AIC: 9075.86 BIC: 9159.275 p: 1 q: 5

## AIC: 9071.572 BIC: 9140.267 p: 2 q: 1

## AIC: 9073.509 BIC: 9147.11 p: 2 q: 2

## AIC: 9075.181 BIC: 9153.689 p: 2 q: 3

## AIC: 9071.24 BIC: 9154.655 p: 2 q: 4

## AIC: 9073.071 BIC: 9161.393 p: 2 q: 5

## AIC: 9073.323 BIC: 9146.924 p: 3 q: 1

## AIC: 9074.454 BIC: 9152.962 p: 3 q: 2

## AIC: 9076.22 BIC: 9159.635 p: 3 q: 3

## AIC: 9072.999 BIC: 9161.32 p: 3 q: 4

## AIC: 9071.952 BIC: 9165.18 p: 3 q: 5

## AIC: 9075.012 BIC: 9153.52 p: 4 q: 1

## AIC: 9076.385 BIC: 9159.8 p: 4 q: 2

## AIC: 9078.449 BIC: 9166.771 p: 4 q: 3

## AIC: 9074.308 BIC: 9167.536 p: 4 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9071.216 BIC: 9169.351 p: 4 q: 5

## AIC: 9076.413 BIC: 9159.828 p: 5 q: 1

## AIC: 9077.747 BIC: 9166.069 p: 5 q: 2

## AIC: 9075.244 BIC: 9168.472 p: 5 q: 3

## AIC: 9076.712 BIC: 9174.847 p: 5 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

## AIC: 9069.659 BIC: 9172.701 p: 5 q: 5

for(i in 1:5) {
  for (j in 1:5){
    simple_model <- arima(diff_relative_humidity_3pm[0:1000], order=c(i,0,j))
    cat()
    cat("AIC:", AIC(simple_model), "BIC:", BIC(simple_model), "p:", i, "q:", j,"\n")
    cat()
  }
}

## AIC: 9074.406 BIC: 9094.037 p: 1 q: 1 
## AIC: 9076.126 BIC: 9100.665 p: 1 q: 2 
## AIC: 9077.978 BIC: 9107.425 p: 1 q: 3 
## AIC: 9079.405 BIC: 9113.759 p: 1 q: 4 
## AIC: 9079.756 BIC: 9119.018 p: 1 q: 5 
## AIC: 9075.475 BIC: 9100.014 p: 2 q: 1 
## AIC: 9077.663 BIC: 9107.109 p: 2 q: 2 
## AIC: 9073.938 BIC: 9108.293 p: 2 q: 3 
## AIC: 9068.856 BIC: 9108.118 p: 2 q: 4 
## AIC: 9069.232 BIC: 9113.402 p: 2 q: 5 
## AIC: 9076.505 BIC: 9105.952 p: 3 q: 1 
## AIC: 9078.539 BIC: 9112.893 p: 3 q: 2 
## AIC: 9074.797 BIC: 9114.059 p: 3 q: 3 
## AIC: 9080.866 BIC: 9125.036 p: 3 q: 4

## Warning in arima(diff_relative_humidity_3pm[0:1000], order = c(i, 0, j)):
## possible convergence problem: optim gave code = 1

## AIC: 9080.726 BIC: 9129.803 p: 3 q: 5

## Warning in log(s2): NaNs produced

## AIC: 9078.304 BIC: 9112.658 p: 4 q: 1 
## AIC: 9079.957 BIC: 9119.219 p: 4 q: 2 
## AIC: 9076.99 BIC: 9121.159 p: 4 q: 3 
## AIC: 9078.994 BIC: 9128.071 p: 4 q: 4 
## AIC: 9078.601 BIC: 9132.587 p: 4 q: 5 
## AIC: 9079.83 BIC: 9119.092 p: 5 q: 1 
## AIC: 9081.592 BIC: 9125.761 p: 5 q: 2 
## AIC: 9078.971 BIC: 9128.049 p: 5 q: 3 
## AIC: 9080.99 BIC: 9134.975 p: 5 q: 4 
## AIC: 9081.844 BIC: 9140.737 p: 5 q: 5

complex_model_lag2 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,], order=c(5, 0, 4))

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1

average_model_lag1 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(5, 0, 5))

## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1

simple_model <- arima(diff_relative_humidity_3pm[0:1000], order=c(2, 0, 4))

AIC(complex_model_lag2)

## [1] 9064.824

BIC(complex_model_lag2)

## [1] 9197.279

AIC(average_model_lag1)

## [1] 9069.659

BIC(average_model_lag1)

## [1] 9172.701

AIC(simple_model)

## [1] 9068.856

BIC(simple_model)

## [1] 9108.118

#Residuals are independent
complex_resids <- residuals(complex_model_lag2)
Box.test(complex_resids,lag=3)

## 
##  Box-Pierce test
## 
## data:  complex_resids
## X-squared = 0.0013766, df = 3, p-value = 1

plot(complex_resids, ylab='Residuals',main='Residual plot for complex model 2')

#Residuals are independent
complex_resids <- residuals(complex_model_lag2)
Box.test(complex_resids,lag=3)

## 
##  Box-Pierce test
## 
## data:  complex_resids
## X-squared = 0.0013766, df = 3, p-value = 1

plot(complex_resids, ylab='Residuals',main='Residual plot for complex model 2')

pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])

ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"),
        ylab='relative_humidity_3pm', main='Actual Vs Predicted')
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)

pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])

ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"),
        ylab='relative_humidity_3pm', main='Actual Vs Predicted')
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)